#######################
#######################
#######################
###Correlation Plots###
#######################
#######################
#######################
library(ggplot2)
library(foreign)
library(zoo)
library(lubridate)

#Set directory
setwd("~/OneDrive - Indiana University/FromGoogle/DiseaseFood/Plots/")

#######
#######
##All##
#######
#######
##Subset data without variables of interest 
plot.dat <- dis.dat[,c("outbreak_event","flus_all","cattle_meat_tonnes1000000","chicken_meat_tonnes1000000","pigs_meat_tonnes1000000","t","year")]


############
###Cattle###
#Correlations
cor(plot.dat$outbreak_event, plot.dat$cattle_meat_tonnes1000000)
#Trendline
lm.pred <- lm(outbreak_event ~ cattle_meat_tonnes1000000, data=plot.dat)
#Confidence intervals
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ cattle_meat_tonnes1000000, data=plot.dat)
se.fit <- sqrt(diag(X %*% V %*% t(X)))
predframe <- with(plot.dat, data.frame(cattle_meat_tonnes1000000,
                                outbreak_event=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

p <- ggplot(plot.dat, aes(x = cattle_meat_tonnes1000000, y = outbreak_event))
#Plot
jpeg("cattleallcor.jpeg", width = 6, height = 6, units = 'in', res = 500)
p +  labs(x="Meat production (tonnes, 1Ms) -- beef", y="Number of outbreaks") +
  geom_point()+ #ylim(162,172) +
  annotate("text", x = 3, y = 250, label = "Correlation = 0.353") +
  geom_line(data=predframe, size=0.5, color="black") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()


#############
###Chicken###
#Correlations
cor(plot.dat$outbreak_event, plot.dat$chicken_meat_tonnes1000000)
#Trendline
lm.pred <- lm(outbreak_event ~ chicken_meat_tonnes1000000, data=plot.dat)
#Confidence intervals
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ chicken_meat_tonnes1000000, data=plot.dat)
se.fit <- sqrt(diag(X %*% V %*% t(X)))
predframe <- with(plot.dat, data.frame(chicken_meat_tonnes1000000,
                                       outbreak_event=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

p <- ggplot(plot.dat, aes(x = chicken_meat_tonnes1000000, y = outbreak_event))
#Plot
jpeg("chickenallcor.jpeg", width = 6, height = 6, units = 'in', res = 500)
p +  labs(x="Meat production (tonnes, 1Ms) -- chicken", y="Number of outbreaks") +
  geom_point()+ #ylim(162,172) +
  annotate("text", x = 7.3, y = 250, label = "Correlation = 0.335") +
  geom_line(data=predframe, size=0.5, color="black") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()



#########
##Pigs###
#Correlations
cor(plot.dat$outbreak_event, plot.dat$pigs_meat_tonnes1000000)
#Trendline
lm.pred <- lm(outbreak_event ~ pigs_meat_tonnes1000000, data=plot.dat)
#Confidence intervals
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ pigs_meat_tonnes1000000, data=plot.dat)
se.fit <- sqrt(diag(X %*% V %*% t(X)))
predframe <- with(plot.dat, data.frame(pigs_meat_tonnes1000000,
                                       outbreak_event=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

p <- ggplot(plot.dat, aes(x = pigs_meat_tonnes1000000, y = outbreak_event))
#Plot
jpeg("pigsallcor.jpeg", width = 6, height = 6, units = 'in', res = 500)
p +  labs(x="Meat production (tonnes, 1Ms) -- pork", y="Number of outbreaks") +
  geom_point()+ #ylim(162,172) +
  annotate("text", x = 28, y = 250, label = "Correlation = 0.346") +
  geom_line(data=predframe, size=0.5, color="black") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()

########
########
##Flus##
########
########

############
###Cattle###
#Correlations
cor(plot.dat$flus_all, plot.dat$cattle_meat_tonnes1000000)
#Trendline
lm.pred <- lm(flus_all ~ cattle_meat_tonnes1000000, data=plot.dat)
#Confidence intervals
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ cattle_meat_tonnes1000000, data=plot.dat)
se.fit <- sqrt(diag(X %*% V %*% t(X)))
predframe <- with(plot.dat, data.frame(cattle_meat_tonnes1000000,
                                       flus_all=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

p <- ggplot(plot.dat, aes(x = cattle_meat_tonnes1000000, y = flus_all))
#Plot
jpeg("cattleflucor.jpeg", width = 6, height = 6, units = 'in', res = 500)
p +  labs(x="Meat production (tonnes, 1Ms) -- beef", y="Number of flu outbreaks") +
  geom_point()+ #ylim(162,172) +
  annotate("text", x = 3, y = 240, label = "Correlation = 0.310") +
  geom_line(data=predframe, size=0.5, color="black") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()


#############
###Chicken###
#Correlations
cor(plot.dat$flus_all, plot.dat$chicken_meat_tonnes1000000)
#Trendline
lm.pred <- lm(flus_all ~ chicken_meat_tonnes1000000, data=plot.dat)
#Confidence intervals
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ chicken_meat_tonnes1000000, data=plot.dat)
se.fit <- sqrt(diag(X %*% V %*% t(X)))
predframe <- with(plot.dat, data.frame(chicken_meat_tonnes1000000,
                                       flus_all=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

p <- ggplot(plot.dat, aes(x = chicken_meat_tonnes1000000, y = flus_all))
#Plot
jpeg("chickenflucor.jpeg", width = 6, height = 6, units = 'in', res = 500)
p +  labs(x="Meat production (tonnes, 1Ms) -- chicken", y="Number of flu outbreaks") +
  geom_point()+ #ylim(162,172) +
  annotate("text", x = 7.3, y = 240, label = "Correlation = 0.293") +
  geom_line(data=predframe, size=0.5, color="black") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()



#########
##Pigs###
#Correlations
cor(plot.dat$flus_all, plot.dat$pigs_meat_tonnes1000000)
#Trendline
lm.pred <- lm(flus_all ~ pigs_meat_tonnes1000000, data=plot.dat)
#Confidence intervals
pp.pred <- predict(lm.pred)
V <- vcov(lm.pred)
X <- model.matrix(~ pigs_meat_tonnes1000000, data=plot.dat)
se.fit <- sqrt(diag(X %*% V %*% t(X)))
predframe <- with(plot.dat, data.frame(pigs_meat_tonnes1000000,
                                       flus_all=pp.pred,lwr=pp.pred-1.96*se.fit,upr=pp.pred+1.96*se.fit))

p <- ggplot(plot.dat, aes(x = pigs_meat_tonnes1000000, y = flus_all))
#Plot
jpeg("pigsflucor.jpeg", width = 6, height = 6, units = 'in', res = 500)
p +  labs(x="Meat production (tonnes, 1Ms) -- pork", y="Number of flu outbreaks") +
  geom_point()+ #ylim(162,172) +
  annotate("text", x = 28, y = 240, label = "Correlation = 0.310") +
  geom_line(data=predframe, size=0.5, color="black") +
  geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)
dev.off()



